Optimal  Control  of  Piezoceramic  Actuators 

Jinghua  Zhong\  Stefan  Seeleckeb,  Ralph  C.  Smithc,  Christof  Btiskensd 


abMechanical  and  Aerospace  Engineering,  NC  State  University,  Raleigh,  NC,  USA  27695 
cCenter  for  Research  in  Scientific  Computation,  NC  State  University,  Raleigh,  NC,  USA  27695 
dDepartment  of  Applied  Mathematics,  University  of  Bayreuth,  Bayreuth,  Germany  95440 


ABSTRACT 

This  paper  presents  the  first  results  of  an  optimal  control  approach  to  piezoceramic  actuators.  A  one-dimensional  free 
energy  model  for  piezoceramics  recently  proposed  by  Smith  and  Seelecke  is  briefly  reviewed  first.  It  is  capable  of 
predicting  the  hysteretic  behavior  along  with  the  frequency-dependence  present  in  these  materials.  The  model  is 
implemented  into  an  optimal  control  package,  and  two  exemplary  cases  are  simulated  to  illustrate  these  features  and  the 
potential  of  the  method. 
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1.  INTRODUCTION 

Actuators  from  piezoceramic  materials  have  established  themselves  in  a  large  number  of  applications  ranging  from 
active  vibration  control  to  nano-scale  positioning  tasks.  This  is  due  to  their  high-frequency  response  behavior  and  their 
essentially  infinite  resolution.  In  particular,  in  the  area  of  nano-scale  positioning,  e.g.,  scanning  microscopy  applications, 
optical  fiber  alignment,  or  general  nano-manufacturing  tasks,  one  is  interested  in  combining  these  two  attractive  features. 
Even  though  tremendous  progress  has  been  made  in  the  area  of  hardware  development  to  enable  higher  speed  in  sensing 
and  measuring,  it  appears  that  the  control  algorithms  represent  the  bottleneck  with  respect  to  improved  performance. 
Most  controllers  are  based  on  conventional  PID  feedback  methods,  which  have  originally  been  developed  for  linear 
systems.  Piezoceramics,  however,  can  be  considered  to  behave  linearly  only  to  about  15%  of  their  maximum  field 
strength.  If  one  is  to  make  more  efficient  use  of  the  actuator,  one  has  to  account  for  the  highly  nonlinear  and  hysteretic 
material  behavior  and  incorporate  it  into  the  control  algorithm. 

A  first  step  in  the  direction  of  improved  control  algorithms  is  a  method  known  as  adaptive  feedforward  control.  It  is 
based  on  the  representation  of  the  system  and  the  control  as  finite  impulse  response  filters  (FIR),  the  coefficients  of 
which  are  determined  by  least  mean  square  algorithms.  A  detailed  overview  of  its  application  to  the  control  of  smart 
structures  is  given  in  [1],  see  [2]  and  [3]  for  a  general  theory  of  adaptive  filters  and  design.  The  method,  however,  is 
restricted  to  linear  systems  or  systems  that  can  be  linearized  and  thus  is  only  applicable  in  the  low  drive  level  range.  A 
next  step  is  the  introduction  of  non-linear  models  into  the  feedforward  control.  In  [4],  a  feedforward  model-reference 
control  method  has  been  applied  to  improve  the  scanning  accuracy  in  a  scanning  tunneling  microscope.  Three  different 
models  have  been  used,  however,  all  of  them  imply  that  the  hysteresis  non-linearity  has  either  local  memory  or 
symmetrical  behavior,  both  of  which  do  not  exactly  represent  the  hysteretic  behavior  of  piezoceramic  actuators.  Also,  a 
pure  feedforward  method  in  a  non-linear  environment  is  always  prone  to  drift  away  from  the  desired  behavior,  so  a 
logical  improvement  is  a  hybrid  method,  which  couples  a  non-linear  model  with  a  feedback  controller.  This  method  is 
also  known  as  inverse  compensator  method,  and  it  has  been  successfully  applied  to  the  control  of  piezoceramic  actuators 
in  [5]  and  [6],  magnetostrictive  transducers  in  [7],  and  shape  memory  alloys  in  [8],  [9]  and  [10].  A  general  overview  can 
be  found  in  [11]  and  [12].  All  these  methods  have  in  common  that  they  try  to  address  robustness  and  stability  issues,  but 
none  of  them  takes  optimality  into  account.  Optimal  control  is  attractive  as  it  aims  to  minimize  some  cost  functional  like 
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energy  consumption  or  final  time.  Some  attempts  toward  this  end  have  been  made  for  shape  memory  alloys  and 
magnetostrictive  materials.  In  [13],  a  SMA  hysteresis  model  is  implemented  into  an  optimal  control  code  to  control  the 
shape  of  an  adaptive  beam  actuated  by  a  SMA  wire.  In  [14]  and  [15],  the  authors  compute  an  optimal  control  for  a 
magnetostrictive  actuator  based  on  an  infinite  time  horizon  optimal  control  problem.  Due  to  the  large  computational 
effort,  these  approaches  have  up  to  only  recently  been  confined  to  offline  simulations.  In  [16],  an  approach  for  a 
magnetostrictive  actuator  is  presented,  which  combines  perturbation  feedback  control  with  a  nonlinear  optimal  control 
method,  thus  allowing  an  online  implementation.  In  a  series  of  papers  [17],  [18],  [19]  and  [20],  a  similar  method  has 
been  developed  for  the  optimal  control  of  shape  memory  alloy  actuators.  It  is  based  on  the  concept  of  a  parametric 
optimal  control  process  combined  with  a  sensitivity  analysis,  and  it  reduces  the  computation  time  for  an  optimal  control 
process  to  the  order  of  nanoseconds.  Thus,  it  has  become  clearly  real-time  capable,  and,  in  particular  through  the 
combination  with  a  strongly  physics-oriented  model  for  the  material  behavior,  the  method  can  be  regarded  as  the  basis 
for  the  development  of  extremely  efficient  control  algorithms  for  active  materials  actuators. 

Due  to  their  close  resemblance  in  microstructure  and  similarity  in  macroscopic  behavior  such  as  hysteresis,  non-linearity 
and  Curie  temperature,  SMAs,  piezoceramics  and  magnetostrictive  materials  can  be  described  in  a  unified  way  [21].  In 
[22],  Smith  and  Seelecke  proposed  a  free-energy-based  model  for  homogeneous  single  crystal  piezoceramics,  which  is  a 
direct  extension  of  the  Muller-Achenbach-Seelecke  theory  for  shape  memory  alloys  [13,18].  The  effects  of  non-uniform 
lattices,  variations  in  effective  fields,  and  polycrystallinity  are  introduced  by  incorporating  appropriate  distributions  in 
the  free  energy  formulations.  The  physical  parameters  in  the  model  greatly  alleviate  the  difficulty  to  accommodate 
changing  operating  conditions  for  traditional  Preisach  models,  which  typically  use  a  large  number  of  unphysical 
parameters  (see  [23]  for  details). 

In  this  paper,  we  combine  the  above  ideas  to  present  an  implementation  of  a  piezoceramic  actuator  model  into  an  optimal 
control  code.  The  implementation  does  not  address  real-time  aspects  yet,  but  as  a  first  step,  focuses  on  the  basic  off-line 
control  problem.  We  briefly  review  the  model  in  section  2,  and  based  on  a  preliminary  implementation  of  the  concept  of 
distributions,  we  illustrate  the  model’s  hysteretic  properties  and  in  particular  its  capability  of  reproducing  the  frequency- 
dependent  behavior  observed  in  these  materials.  In  the  final  section,  we  discuss  the  implementation  into  an  optimal 
control  code  and  present  two  non-trivial  control  problems  addressing  the  material  properties  mentioned  above. 

2.  A  MODEL  FOR  PIEZOCERAMIC  MATERIALS 


2.1  Single  crystal 

This  section  introduces  the  model,  which  is  used  in  the  subsequent  optimal  control  problems.  Numerous  models  have 
been  proposed  in  the  past  decades  to  model  hysteresis,  most  of  which  are  motivated  by  research  in  magnetic  materials. 
For  piezoceramics,  microscopic  models  at  the  grain  or  lattice  level  are  rarely  employed  for  control  purposes  due  to  the 
large  number  of  parameters  and  states  involved.  Macroscopic  models  are  typically  low-order  and  based  on  empirical 
principles,  e.g.  the  Preisach  model  [26]  and  its  variations  [27,28],  but  the  absence  of  physical  parameters  makes  it 
difficult  to  incorporate  frequency,  temperature  or  load  dependencies  and  adapt  to  changing  operating  conditions.  Also, 
most  of  these  models  only  consider  part  of  the  hysteresis  loop  with  a  unipolar  field  [6],  which  limits  the  available  stroke 
of  the  actuators.  Semi-macroscopic  (mesoscopic)  models,  including  the  model  used  in  this  paper,  are  derived  from 
energy  principles  and  try  to  combine  microscopic  mechanisms  with  macroscopic  averages  to  obtain  physical  parameters 
[29,30].  They  usually  provide  a  good  compromise  between  capturing  the  physical  properties  of  the  material  and  being 
computationally  efficient,  so  that  they  represent  a  good  choice  for  the  implementation  into  a  control  algorithm. 

In  this  paper,  we  use  the  one-dimensional  energy  model  proposed  by  Smith  and  Seelecke  for  piezoceramics.  It  is 
developed  with  the  assumption  that  dipoles  in  a  single  crystal  with  uniform  lattice  spacing  have  two  orientations,  which 

are  denoted  by  (+)  and  (— )  ,  depending  on  whether  the  77 41  ion  is  displaced  upward  or  downward  with  respect  to  the 
center  of  symmetry  (see  Figure  1 ). 

For  a  typical  mesoscopic  domain  consisting  of  a  number  of  equally  poled  lattice  cells,  the  model  assumes  the  following 
multiparabolic  representation  for  the  Helmholtz  free  energy  density  If/  as  a  function  of  the  polarization  P : 
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HP)  = 


(i) 


\EX{P  +  Pr f  ,P<-PX 
<  \Ex(P-Pt)2  ,P  >  P\  . 

|£1(P1-JPr)[^-JPr],|JP|<JP1 

V  lii 

In  Eq.  ( 1 ),  PT  ,  P{  and  Ex  are  related  to  physical  properties,  which  can  be  extracted  from  an  electric  field  vs. 

polarization  diagram.  PT  denotes  the  remnant  polarization,  Px  the  onset  of  the  switching  process,  and  Ex  is  the 
reciprocal  of  the  slope  of  the  ( E ,  P)  curve  (see  the  bottom  row  of  Figure  2). 

The  switching  processes  in  the  material  can  be  interpreted  as  phase  transformations  with  the  polarization  as  the  order 
parameter.  The  kinetics  of  these  phase  transformations  are  determined  by  the  electric-field-dependent  barriers  and 
minima  in  the  Gibbs  free  energy,  which  is  related  to  the  Helmholtz  free  energy  in  Eq.  (1)  through  the  Legendre 
transformation  G  =  lj/  —  EP  .  This  dependence  is  illustrated  in  Figure  2  for  increasing  field  values  (from  left  to  right) 
together  with  the  corresponding  polarization  vs.  electric-field  behavior. 


Figure  1  PbTiCfi  crystal  lattice  cells  with  positive  and  negative  polarization 


V(P.T)  =  G(0,P,T) 


(b) 

Figure  2  (a)  Helmholtz  energy  i //  and  Gibbs  energy  G  for  increasing  field  E\ 
(b)  Polarization  P  as  a  function  of  E  for  a  single  crystal  with  uniform  lattice. 
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Due  to  the  ever-present  thermal  activation,  all  domains  fluctuate  about  their  equilibrium  polarization  values.  From 
Boltzmann  principles,  the  probability  of  attaining  an  energy  level  G  is 

ju(G)  =  Ce~GlkT ,  (2) 

where  k  denotes  Boltzmann's  constant.  The  expected  average  polarizations  due  to  domains  with  positive  and  negative 
polarizations  are  then  given  by 


(P+)  = 
(P-) 


r  pe-G<E,r;nvnikrdp 
_ 

£ 

r> 


-G(E,P,T)VD!kT 


,  and 


dP 


-G(E,P,T)Vd  /  kT 


-G(E,P,T)VDlkT 


dP 


and  the  average  polarization  for  the  whole  actuator  is 

P  =  xJP,)  +  x_(P_), 


(3) 


(4) 


(5) 


where  X+  and  X_  denote  the  fractions  of  domains  with  positive  or  negative  polarization  respectively.  These  phase 
fractions  evolve  in  time  according  to  the  differential  equations 


i+  =  ~P+-X+  +  p_+x_ 

(6) 

X-  =  ~P-+X_  +  P+_X+ 

which  can  be  simplified  to 

*+  =  ~P+- X+  +  p_+  (1  -  X+  )  (7) 


through  the  identity  X+  +  X_  =  1 .  Here  p+  and  p_+  denote  the  transition  probabilities  of  a  domain  switching  from 
one  polarization  to  the  other.  They  are  computed  by 


1 

**  r e-G^px)vDikTdP 

g-G^E-P^TXTWolkT 
Tx  r''  e-G(E,p,r}vnlkT dp 


(8) 


where  Tx  is  the  relaxation  time  of  the  material  and  V D  denotes  the  activation  volume  of  the  domain.  The  ratio  V D  /  kT 

determines  the  fluctuation  of  polarization  about  the  equilibrium  values  due  to  thermal  activation,  and  together  with  T  , 

this  ratio  dictates  the  frequency  dependence  of  the  hysteresis  loops.  The  solution  of  equation  (5)  specifies  the  hysteretic 
relation  between  E  and  P,  as  is  depicted  in  Figure  2. 
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2.2  Polycrystalline  materials 


As  the  above  model  describes  perfect  single  crystal  behavior,  realistic  actuators  are  neither  perfect  nor  are  they  single 
crystals.  Rather  they  exhibit  inhomogeneities  due  to  lattice  imperfections,  impurities,  grain  boundaries,  etc.,  which  lead 
to  variations  in  the  energy  barriers,  depending  on  the  location  of  the  domain  under  consideration.  Following  again  the 
ideas  outlined  in  [22]  and  [24],  we  take  these  variations  into  account  by  introducing  a  distribution  in  the  energy  barriers. 

In  particular,  we  consider  the  lattice  parameter  Ex  =  Ex  (  PT  —  Px  )  to  be  normally  distributed  with  mean  E  .  The  total 
polarization  is  then  given  by 


P(E)  =[P(E,  Ex)f(Ex)dEv  (9) 

with  the  density 

f(Ex)  =  Ce~(S'~W'f,b .  (10) 

Other  aspects  such  as  variations  in  the  effective  field  can  be  addressed  similarly  by  introducing  additional  distributions 
for  the  corresponding  parameters  [22].  To  simplify  the  algorithm  and  reduce  computational  cost  for  the  solution  of  the 
optimal  control  problem,  only  the  distribution  of  energy  barriers  is  considered  and  implemented  in  this  paper. 

2.3  Numerical  implementation 

In  [18],  the  integral  in  Eq.  (9)  is  numerically  solved  by  a  Gaussian  quadrature  method.  A  set  of  single  crystal  ODEs  must 
be  solved  at  each  of  the  abscissae  points,  so  that  the  number  of  Gauss  points  determines  the  size  of  the  system.  Obtaining 
a  smooth  hysteresis  curve  thus  requires  the  solution  of  a  large  number  of  ODEs,  which  severely  reduces  the 
computational  speed.  In  [24],  the  authors  suggest  an  efficient  scheme  to  calculate  the  overall  polarization  from  the 
switching  of  individual  domains  by  utilizing  algebraic  matrix  operations,  but  it  is  valid  only  for  the  limiting  case  of  zero 
thermal  activation,  and  is  not  able  to  account  for  effects  related  to  different  loading  rates.  Within  the  context  of  time- 
dependent  processes,  we  are  particularly  concerned  with  the  model's  capability  of  reproducing  the  frequency  dependent 
response  correctly.  In  order  to  retain  this  capability,  we  choose  yet  an  alternative  way  to  implement  the  distributions, 
which  at  the  same  time,  proves  to  be  highly  computationally  efficient. 

To  reduce  computational  cost,  we  implement  the  distribution  into  the  ODE  by  parameterizing  the  process  of  the  phase 
transition.  Instead  of  calculating  many  domains  simultaneously,  we  consider  a  representative  domain,  which  is  most 

likely  to  transform  next  at  a  given  value  of  the  phase  fraction  and  the  prescribed  electric  field.  The  lattice  parameter  E x 
for  this  domain  is  calculated  from  the  current  phase  fraction  and  loading  history.  This  way,  only  one  ODE  needs  to  be 
solved.  A  preliminary  implementation  of  this  concept,  which  is  used  for  the  optimal  control  code,  correctly  predicts  the 
outer  hysteresis  loop,  although  further  development  is  necessary  to  retain  the  correct  inner  loops.  Details  will  be  given  in 
a  forthcoming  paper  [32].  Figure  3  illustrates  the  model's  capability  of  capturing  the  frequency  dependence  for  two 
typical  values.  Note  the  dramatic  impact  on  the  shape  of  outer  and  inner  hysteresis  loops.  The  parameter  values  used  for 

the  simulations  throughout  the  paper  are:  T  —  273.0  K,  VD  =  1.0x10  22m3,  Tx  =  1.0x10  2s_1,  Px  =0.275C/m2, 
PT  =0.26  C/m2,  Ex  =  6.67  X 107  V-m/C,  Z7  =  2.0xl012,  and  E  =1.0xl06V/m. 

3.  OPTIMAL  CONTROL 


3.1  Optimal  control  problem 

In  the  following  section,  we  will  study  the  solution  of  the  optimal  control  problems.  For  simplicity,  we  consider  a  one¬ 
dimensional  piezoceramic  stack  actuator  in  the  absence  of  mechanical  loads,  and.  instead  of  tracking  the  displacement, 
we  consider  the  polarization  as  the  quantity  that  can  be  prescribed. 


5 


Thus  we  prescribe  the  set-point  polarization  Pet(t )  ,  and  solve  for  the  electric  field  function  E{t)  ,  which  minimizes  the 
following  cost  functional 


J=fo[P(t)-Pset(t)fdt. 


(11) 


PZT  Simulation  -  Electric  Field 


PZT  Simulation  -  P  vs  E 


Electric  Field  (V/m) 


PZT  Simulation  -  Electric  Field 


PZT  Simulation  -  P  vs  E 


Electric  Field  (V/m) 


Figure  3  Left:  hysteresis  curve  at  0.125Flz,  Right:  hysteresis  curve  at  125FIz 


This  ensures  that  the  set  point  is  reached  in  the  shortest  time  possible  and  subsequently  maintained.  The  minimization  is 
subject  to  the  following  constraints 


X  =  f{x(t),E(t)l 


x(0)  =  x0, 

\m\  <  Em 


and 


(12) 


Equation  (12)i  represents  the  model  equation,  which  describes  the  switching  processes  in  the  actuator,  Eq.  (12)2  are  the 
initial  conditions  of  the  process,  and  Eq.  ( 12)3  are  box  constraints  for  the  control,  representing  the  fact  that  only  a  finite 
electric  field  can  be  applied  to  the  actuator. 
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For  the  numerical  solution,  we  have  implemented  the  model  from  the  previous  section  into  the  optimal  control  package 
NUDOCCCS  developed  by  Biiskens  [25],  Based  on  a  direct  approach,  NUDOCCCS  discretizes  the  control  and 
transforms  the  original  control  problem  into  a  nonlinear  optimization  problem.  The  control  function  (the  electric  field 
E(t ) )  can  be  either  assumed  piecewise  linear  or  interpolated  by  higher  order  splines.  We  have  used  the  piecewise  linear 
option  for  simplicity  and  faster  calculation. 

The  differential  equations  for  the  state  variables  are  integrated  by  an  implicit  Runga-Kutta  scheme  (RADAU  Ha)  using 
the  efficient  and  robust  RADAU5  routine  from  the  book  of  Hairer  and  Wanner  [31].  NUDOCCCS  then  solves  the 
resulting  NLP  problem  with  a  sequential  quadratic  programming  method. 

3.2  Case  1:  Single-level  set  point 

We  start  with  a  simple  case,  which  allows  illustrating  the  rate-dependent  behavior  of  the  actuator.  From  an  initial 
polarization  P(0)  =  0  C/m2,  the  set  point  rises  linearly  to  a  new  level  Pset  =  0.2  C/m2  at  a  rate  of  500  C/m2/sec  within 

4  ms.  The  constraint  E  has  been  set  to  6  MV/m. 


The  four  plots  in  Figure  4  show  the  time  evolutions  of  the  phase  fraction  x(t)  and  polarization  Pi  t)  along  with  the 
value  of  the  cost  functional  and  the  optimal  electric  field  E(t)  (from  upper  left  to  lower  right).  Figure  5  shows  the 
hysteretic  behavior  of  the  actuator  in  a  P-E  plot. 


Time  (sec) 


Time  (sec) 


Figure  4  Optimal  control  of  piezoceramic  polarization  (Case  1) 
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The  polarization  calculated  from  the  optimal  solution  follows  the  set-point  function  closely,  especially  in  the  middle  of 
the  rising  process,  which  can  be  seen  from  the  cost  functional  increasing  only  very  little  during  this  phase  of  the  process. 
Most  of  the  error  occurs  near  the  beginning  and  the  end  of  the  transition.  In  the  mean  time,  the  phase  fraction  increases 
monotonously  from  0.5  to  0.83.  The  electric  field,  however,  exhibits  some  interesting  features  due  to  the  nonlinear  and 
hysteretic  behavior  of  piezoceramics  at  different  loading  rates. 

When  a  positive  electric  field  is  applied  to  the  piezoceramic  actuator,  the  increase  in  polarization  contains  two  different 
terms.  The  first  term  is  proportional  to  the  applied  field  and  time  independent.  The  second  term  is  governed  by  the 
evolution  of  the  phase  fraction  and  as  such  is  time  dependent.  The  model  automatically  accounts  for  this  behavior 
through  a  suitable  choice  of  relaxation  time  and  activation  volume,  and  consequently,  the  piezoceramic  actuator  exhibits 
different  shapes  of  hysteresis  at  different  loading  frequencies. 

As  we  can  see  from  Figure  4,  the  electric  field  starts  with  a  sharp  increase,  which  is  equivalent  to  a  loading  frequency  of 
1250  Hz,  then  continues  to  increase  at  a  roughly  10  times  slower  rate,  until  it  reaches  the  peak  of  3.3  MV/m.  At  this 
electric  field  value,  the  polarization  has  almost  reached  the  required  set  point,  but  the  phase  transition  has  not  yet 
completed.  Some  domains  of  negative  phase  continue  to  change  into  positive  phase,  which  will  increase  the  overall 
polarization  of  the  actuator.  To  compensate  for  this  increase,  the  electric  field  then  decreases  asymptotically  to  about  1.9 
MV/m.  This  process  is  better  illustrated  in  Figure  5,  which  shows  the  optimal  solution  along  with  hysteresis  loops  at 
three  different  loading  frequencies. 


Figure  5  P-E  behavior  of  the  optimal  solution  (Case  1) 

In  Figure  5,  the  actuator  first  follows  the  hysteresis  curve  at  the  highest  frequency  (the  widest  loop),  then  as  the  rate  of 
the  electric  field  decreases,  it  overlaps  with  the  curve  at  a  lower  loading  frequency,  and  finally  relaxes  to  a  point  on  the 
quasistatic  outer  loop  (0.125  Hz).  At  a  frequency  of  0.125  Hz,  the  effect  from  the  relaxation  time  of  the  domains  is 
negligible  because  the  domains  now  have  enough  time  to  complete  the  phase  transition  as  the  field  increases.  If  no 
optimal  control  is  applied,  the  same  polarization  set  point  can  be  achieved  by  slowly  increasing  the  electric  field.  The 
actuator  would  then  follow  the  hysteresis  curve  of  0.125  Hz. 

3.3  A  more  complex  set-point  function 

Due  to  hysteresis,  the  same  polarization  may  require  a  different  electric  field  if  the  loading  history  is  different.  We 
proceed  with  the  prescription  of  a  more  complex  set-point  function  to  verify  the  model's  ability  to  predict  the  polarity  of 
the  electric  field  needed  for  a  positive  polarization  with  different  history.  The  set-point  polarization  first  increases  from 


0.0  C/m2  to  0.1  C/m2  linearly  at  a  rate  of  500  C/m2/sec,  holds  for  a  period  of  time,  increases  again  to  0.2  C/m2  at  the 
same  rate,  holds  the  value,  and  then  decreases  back  to  0. 1  C/m2  at  a  rate  of  -500  C/m2/sec.  In  the  simulation,  the  control 

process  has  a  total  duration  of  4  ms,  and  is  discretized  using  81  points  (NDISKRET=81  in  NUDOCCCS).  E  is  set  to 
6  MV/m.  The  four  plots  in  Figure  6  show  the  time  evolutions  of  the  phase  fraction  x(t)  and  polarization  Pit )  along 
with  the  value  of  the  cost  functional  and  the  necessary  electric  field  E(t )  (from  upper  left  to  lower  right).  Figure  7 
shows  the  hysteretic  behavior  of  the  actuator  in  a  P-E  plot. 


Time  (sec)  Time  (sec) 


Time  (sec)  Time  (sec) 

Figure  6  Optimal  control  of  piezoceramic  polarization  (Case  2) 


In  this  case,  not  surprisingly,  the  model  predicts  similar  shapes  of  optimal  electric  fields  for  the  actuator.  The  electric 
field  increases  first  to  a  high  value,  and  then  drops  back  to  a  lower  value  as  the  piezoceramic  material  relaxes,  forcing 
the  polarization  to  follow  the  set  point.  The  maximum  electric  field  needed  for  the  second  set  point  of  0.2  C/m2  is  almost 
the  same  as  that  of  the  previous  section,  and  the  field  also  decreases  at  the  same  rate  to  the  same  final  value  to  maintain 
the  polarization.  As  Figure  7  shows,  in  following  the  three  set  points,  the  actuator  first  follows  the  hysteresis  loop  of  a 
higher  frequency  but  always  comes  back  to  a  point  on  a  loop  of  low  frequency. 

For  the  first  and  third  set  point,  although  they  have  the  same  polarization  of  0. 1  C/m2,  the  model  predicts  completely 
different  phase  fractions  and  electric  fields.  Obviously,  this  is  due  to  the  difference  in  the  history.  For  the  first  set  point,  a 
positive  electric  field  is  applied,  because  the  initial  positive  phase  fraction  needs  to  be  increased.  To  reach  the  third  set 
point,  the  hysteresis  has  to  be  traversed,  and  a  negative  electric  field  has  to  be  applied  in  order  to  reach  the  same 
polarization  value  as  for  the  first  set  point. 


9 


Figure  7  P-E  behavior  of  the  optimal  solution  (Case  2) 


4.  CONCLUSION 

The  paper  introduced  an  implementation  of  an  energy-based  model  for  the  simulation  of  polycrystalline  piezoceramic 
materials.  The  model  reproduces  the  outer  hysteresis  loops  typically  observed  in  these  materials  correctly,  although 
further  development  is  needed  for  a  better  simulation  of  minor  loops.  It  was  also  shown,  that  the  effect  of  varying 
loading  rate  is  captured,  and  both  these  effects  were  illustrated  to  play  an  important  role  in  typical  control  problems.  The 
model  has  been  implemented  into  an  optimal  control  code,  which  does  not  only  compensate  for  the  non-linear  and 
hysteretic  material  behavior,  but  allows  optimality  criteria  like  speed  of  adjustment  to  be  taken  into  account.  This  first 
implementation  demonstrated  the  feasibility  of  the  concept  by  successfully  computing  offline  solutions,  while  future 
work  will  focus  on  implementing  real-time  concepts  recently  developed  by  the  authors  for  shape  memory  alloys  [16,17]. 
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